Bayesian Empirical Generator

In this section we will show how to utilise the Bayesian empirical generator to estimate the uncertainty of the entries in the empirical generator from finite data. We we will do so by generating a Markov chain from a generator matrix. First let's load a few packages and functions that we've become familiar with over the last few sections.

using Random, MarkovChainHammer, Statistics
using MarkovChainHammer.TransitionMatrix: generator
using MarkovChainHammer.Trajectory: generate
Random.seed!(1234)

Random.TaskLocalRNG()

Observe the use of the Random package to set the seed for the random number generator. This is important for reproducibility of the results. A seed should only be defined once at the beginning of a script.

We will now create a generator matrix Q and generate a markov process.

Q = [-1.0 4/3 2; 1/4 -2.0 1; 3/4 2/3 -3.0]
dt = 0.01
markov_chain = generate(Q, 10000; dt = dt)'

1×10000 adjoint(::Vector{Int64}) with eltype Int64:
 1  1  1  1  1  1  1  1  1  1  1  1  1  …  1  1  1  1  1  1  1  1  1  1  1  1

We first construct the empirical generator as before to serve as a comparison point for the Bayesian empirical generator which will come later.

Qempirical = generator(markov_chain; dt = dt)

3×3 Matrix{Float64}:
 -1.13282    1.10375    2.28898
  0.266546  -1.65563    0.66762
  0.866273   0.551876  -2.9566

We are now ready to introduce the Bayesian Generator matrix. We first load in the BayesianMatrix submodule.

using MarkovChainHammer.BayesianMatrix

This submodule naturally exports a number of structs and functions. We can see the names of these functions by typing the following:

names(MarkovChainHammer.BayesianMatrix)

8-element Vector{Symbol}:
 :BayesianGenerator
 :BayesianMatrix
 :GeneratorParameterDistributions
 :eigen_distribution
 :eigvals_distribution
 :params
 :uninformative_prior
 :unpack

In particular a BayesianGenerator object is exported. We use the BayesianGenerator just like a normal generator in order to construct a BayesianGenerator

Q_bayes = BayesianGenerator(markov_chain; dt = dt)

BayesianGenerator 
prior variance  
3×3 Matrix{Float64}:
 7.03687e13  3.51844e13  3.51844e13
 3.51844e13  7.03687e13  3.51844e13
 3.51844e13  3.51844e13  7.03687e13prior mean  
3×3 Matrix{Float64}:
 -1.0   0.5   0.5
  0.5  -1.0   0.5
  0.5   0.5  -1.0posterior variance  
3×3 Matrix{Float64}:
 0.0185982   0.0609135  0.109155
 0.00442456  0.0913703  0.0318369
 0.0142707   0.0304568  0.140992posterior mean  
3×3 Matrix{Float64}:
 -1.13282    1.10375    2.28898
  0.266546  -1.65563    0.66762
  0.866273   0.551876  -2.9566

This is no longer a regular matrix, but rather a random matrix whose entries are given by a probability distribution consistent with finite sampling from a markov process. In the present context (this will change later when we discuss prior distributions) the mean of the Bayesian matrix reproduces the same matrix as the empirically obtained matrix, as we can check

mean(Q_bayes) - Qempirical

3×3 Matrix{Float64}:
 -2.22045e-16  -2.22045e-16  -1.33227e-15
  2.22045e-16   4.44089e-16   1.11022e-16
  1.11022e-16   2.22045e-16   1.33227e-15

However, we have more than just a mean value for the Bayesian matrix, we also have variances

var(Q_bayes)

3×3 Matrix{Float64}:
 0.0185982   0.0609135  0.109155
 0.00442456  0.0913703  0.0318369
 0.0142707   0.0304568  0.140992

and standard deviations

std(Q_bayes)

3×3 Matrix{Float64}:
 0.136375   0.246807  0.330386
 0.0665173  0.302275  0.178429
 0.11946    0.174519  0.375489

We can check that the mean falls within two standard deviations of the true empirical estimate with high probability

abs.(mean(Q_bayes) - Q) .<  2 * std(Q_bayes)

3×3 BitMatrix:
 1  1  1
 1  1  1
 1  1  1